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ABSTRACT 

The equation describing the secular diffusion of a self-gravitating collisionless system 
induced by an exterior perturbation is derived while assuming that the timescale 
corresponding to secular evolution is much larger than that corresponding to the 
natural frequencies of the system. Its two dimensional formulation for a tepid galactic 
disc is also derived using the epicyclic approximation. Its WKB limit is found while 
assuming that only tightly wound transient spirals are sustained by the disc. It yields 
a simple quadrature for the diffusion coefficients which provides a straightforward 
understanding of the loci of maximal diffusion within the disc. 
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1 INTRODUCTION 

Understanding the secular dynamical evolution of galaxies 
over cosmic time has been a long standing subject of in¬ 
terest. Indeed, self gravitating collisionless systems such as 
galaxies may, over cosmic times, change their kinematical 
structure as they respond secularly to their evolving envi¬ 
ronment, in a manner which depends both on their internal 
orbital structure, but also on how this structure resonates 
with its environment or with itself. It is therefore critical to 
distinguish in the physical properties of galaxies the contri¬ 
butions from the cosmic environment (nurture) and its in¬ 
duced pertubations from the ones coming from the intrinsic 
properties of the galaxies (nature). For thermodynamically 
improbable cold systems such as galactic discs, their gravi¬ 
tational susceptibility should also play a specific role which 
must be taken into account when studying their long term 
evolution. 

To tackle this question, one can rely on numeri¬ 
cal N-body simulations of higher resolutions to take into 
account non-linear physical processes (e.g. Dubois et al. 
2014) or perform idealized well-crafted numerical exper¬ 
iments (Sellwood & Athanassoula 1986; Earn & Sellwood 
1995; Sellwood 2012). With such statistical investigations, 
one can assess the importance and the role of the or¬ 
bital structure of a galactic disc to drive its secular evo¬ 
lution. Angle-action variables (Born 1960; Goldstein 1950; 
Binney & Tremaine 2008) and the matrix method (Kalnajs 
1976) also allow us to take into account the self-gravitating 


amplification of such collisionless systems. From this analyti¬ 
cal framework, one should therefore be able to derive flexible 
qualitative and quantitative equations describing the secular 
dynamics of discs, without relying on the implementation of 
the corresponding demanding numerical models. 

This topic of secular evolution has been addressed via 
the dressed Fokker-Planck equation, where the source of sec¬ 
ular evolution for a self-gravitating system is taken to be po¬ 
tential fluctuations from an external bath, e.g. correspond¬ 
ing to the cosmic environment. Binney & Lacey (1988) com¬ 
puted the first and second-order diffusion coefficients de¬ 
scribing the orbits deviation induced by fluctuations in the 
gravitational potential. Weinberg (1993) showed the impor¬ 
tance of self-gravity on the nonlocal and collective relaxation 
of stellar systems. Weinberg (2001a) and Weinberg (2001b) 
considered the dressed gravitational amplification of Poisson 
shot noise in stellar systems and the impact of the proper¬ 
ties of the noise processes. Ma & Bertschinger (2004) used 
a quasi-linear approach to investigate dark matter diffu¬ 
sion induced by cosmological fluctuations. Pichon & Aubert 
(2006) sketched a time-decoupling approach to solve the col¬ 
lisionless Boltzmann equation and applied it to the statis¬ 
tical study of dynamical flows through dark matter haloes. 
Chavanis (2012b) considered the evolution of homogeneous 
collisionless systems forced by an external perturbation. 
Nardini et al. (2012) also considered the evolution of such 
long-range interacting systems when driven by external 
stochastic forces. 

Using an argument based on timescale decoupling in- 


2 Jean-Baptiste Fouvry, Christophe Pichon, Simon Prunet 


spired from Pichon & Aubert (2006), we present here a care- 
fnl and detailed derivation of the secular resonant dressed 
orbital diffusion equation for a general collisionless self- 
gravitating system undergoing external perturbations ^. We 
then develop this formalism for the secnlar evolution of 
an infinitely thin galactic disc. In order to circnmvent the 
complex direct or analytical calculation of the modes of a 
galactic disc carried only for a small number of disc mod¬ 
els (Zang 1976; Kalnajs 1977; Goodman 1988; Weinberg 
1991; Vauterin & Dejonghe 1996; Pichon & Cannon 1997; 
Evans & Read 1998; Jalali & Hunter 2005), we then rely 
on the WKB approximation (Liouville 1837; Toomre 1964; 
Kalnajs 1965; Lin & Shu 1966) to obtain a tractable alge¬ 
braic expression for both the gravitational snsceptibility of 
the system and the associated diffnsion coefficients. Within 
the realm of this approximation, which shonld apply to 
cold enongh discs, the diffnsion coefficient reduce to simple 
quadratures. 

The paper is organized as follows. Section 2 presents 
a derivation of the general dressed Fokker-Planck equation 
for a pertnrbed self-gravitating collisionless system. Ap¬ 
pendix A provides a complementary derivation based on 
Hamilton’s equation, extending Binney & Lacey (1988) to 
self gravitating systems. Section 3 focuses on razor thin ax- 
isymmetric galactic discs within the WKB approximation. 
Some of the nnderlying calculations are postponed to Ap¬ 
pendixes B and C. Finally, section 4 wraps up. 


2 SECULAR DIFFUSION EQUATION 

The secular diffusion equation aims at describing the long¬ 
term aperiodic evolution of a self-gravitating collisionless 
system, perturbed by exterior potential fluctuations. A typ¬ 
ical application for this formalism is the study of a galactic 
disc undergoing (cosmic) perturbations from its surround¬ 
ing dark matter halo or the secular diffusion of accretion 
streams within the Galactic halo. We will suppose that the 
background gravitational potential of the system is station¬ 
ary and integrable, so that we may always remap the usual 
phase-space coordinates (x,v) to the angle-action coordi¬ 
nates (0, J). This is a strong assumption, as one could imag¬ 
ine situations where the secular evolution breaks symme¬ 
try warranting integrability. The angles 0 are 27r—periodic, 
whereas the actions J are conserved for a few dynamical 
times and secularly drift with cosmic time. 


2.1 Evolution equations 

We consider a stationary Hamiltonian Hq{J), associated to 
a stationary background gravitational potential tpo- We also 
consider a quasi-stationary distribution function Fo{J,t), 
which, at fixed secular time, only depends on the actions 
thanks to Jeans theorem (Jeans 1915). Finally, we suppose 
that an exterior source is perturbing this stationary system, 
so that we can expand the distribution function and the 


^ We also recover a missing 1/2 factor absent from this previous 
work, due to an error in temporal integration bounds. 


Hamiltonian of the system as 


F{J,0,t)^FoiJ,t) + fiJ,0,t), 

H{J, 0, t) = Ho (J) -b t/>'=(J, 0, t) + (J, 0, t), 


where / is the perturbation of the distribution function, t/)® 
is the perturbing exterior potential generated by the ex¬ 
terior source, and is the self-response from the system 
induced by its self-gravity. This decomposition now involves 
two main temporal scales. The shortest scale is the fluc¬ 
tuation timescale, during which Fo{J) may be considered 
constant. The longest timescale corresponds to the secular 
evolution timescale. The perturbations are supposed to be 
small so that / <C Fb and ,%l}‘ ifo. The evolution of the 
collisionless system is then driven by Boltzmann collisionless 
equation which reads 


AF 

dt 


dt 


+ {H,F} = 0, 


( 2 ) 


where {H,F} is the Poisson bracket. Injecting the decom¬ 
position from equation (1) in Boltzmann’s equation (2), we 
obtain 


0 = 


dFo df 
dt dt 




~d0^^ 

dxf^ d^ 
dJ ^ dJ 


dFo 

dJ 

di 

d0 ’ 


dip^ dip" 

( 3 ) 


dJ 


where we have defined the frequencies of the motion on the 
action-torii as 


0 = n = 


dHo 

dJ ■ 


( 4 ) 


In order to derive the corresponding secular equation, we 
perform an angle-average on 0 of equation (3). All terms 
involving a single derivation d/d0 are equal to 0, since the 
angles 0 are 27r-periodic. Moreover, we have Jgd0 f = 0, 
because all the variations independant of 0 are included in 
Fo{J, t). As Fo is independent of 0, we obtain 


dFo 

1 

1 d0 

'd-ip^ dip”' 

df 

dt 

(27r)‘* J 


dJ 


1 

! d0 

'dip” dip”' 

df 


(27r)‘* J 

dJ^dJ 

d0 ’ 


where d is the dimension of the physical space. Using 
Schwartz theorem, this secular diffusion equation can be 
written under the shorter form 

dFo _ _L_^ 

dt dJ 


d0f 


dlip'^ d-pj" 


d0 


( 6 ) 


Equation (6), written as the divergence of a flux, empha¬ 
sizes the fact that during orbital diffusion the total num¬ 
ber of stars is strictly conserved. Recalling that / <C Eq and 
<C 'i/’o, the secular evolution equation (6) shows that 
dFo/dt is in fact a second order term. 

Gorrespondingly, keeping only hrst order-terms in (3), 
we obtain the second diffusion equation for the short 
timescale, which reads 


^ o 9[r+r] 

dt d0 dJ ' d0 


( 7 ) 


This equation describes the evolution of the perturbation 
distribution function / on the fast fluctuating timescale. On 
such timescales, dFo/dJ will be considered as independant 
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of t. The next step is to study the fast fluctuating equa¬ 
tion (7), whose solutions will allow us to estimate the dif¬ 
fusion coefficients for the secular evolution given by equa¬ 
tion (6) and describe the diffusion of the quasi-stationary 
distribution function Fq in action-space. 


2.2 Fourier expansion 


One of the many assets of the angle-action variables is that 
the angles 0 are 27r—periodic allowing us to perform discrete 
Fourier expansions with respect to these variables. We define 
the Fourier transform in angles of a function X{G, J) as 

'X{0,J)= 

yd0X(e, 


Thanks to this transformation, the evolution equation (7) 
takes the form 


dfm 

dt 


-I- im-Fi fm 


— im- 


dFo 

dJ 


[i’L+i’L] = 0. 


(9) 


At this stage, we introduce the assumption of timescale de¬ 
coupling and push the secular time to inhnity. As a conse¬ 
quence, in the upcoming calculations, we will suppose that 
dFo/dJ = cst. Forgetting transient terms and bringing the 
initial time to —oo, to focus only on the forced regime, the 
equation (9) can be solved explicitly, leading to 


dre . (10) 


We define the temporal Fourier transform as 

r+oo 


/H= / 

J —oo 

-1 ^ + oo ^ 

/(i) = ^ / daj/(w)e" 


( 11 ) 


Taking the Fourier transform of equation (9) at the fre¬ 
quency uj, one can write 

= [ fc ( J ,^)+?^( J ,^)1 . ( 12 ) 

uj — mll L J 


2.3 Matrix method 

An important property of this self gravitating system is that 
the perturbed distribution function / is consistent with the 
self-gravitating potential and its associated density p®, 
so that we have 

p‘{x,t)= f dv f{x,v,t) f dv fm{ J,t)e'^'*’. (13) 

•' m •' 

In order to simplify further equation (13), we follow Kalnajs 
matrix method (Kalnajs 1976) and introduce a complete 
biorthonormal basis of potentials and densities and 

p^^\x), such that 

( = 4-kGp^^'> , 

^Jdx[tp^^\x)]*p^'^\x) = —Sp. ^ ^ 


On such a basis, we can write the exterior and the self po¬ 
tentials as 


= ^ap(t) , 

pGN 


pGN 


(15) 


The linearity of Poisson’s equation ensures that this decom¬ 
position translates into p‘‘{x,t) — '^^ap{t)p^^\x). Using the 
biorthogonality of the basis, we multiply equation (13) by 
(a;)]* and integrate over all positions to obtain 

ap{t) = Jdxdv [ip^^\x)]* . (16) 


2.4 Response matrix and self-consistency 

As the transformation {x, v) i—>■ (0, J) is canonical, we have 
da:du = d0dj. The integration on 0 in equation (16) is 
straigthforward since only (a;)]*e'”^'® depends on it, so 
that it becomes 

ap(t) = -( 2 nfY [V’^’(J)]‘ • (17) 

m 

Using the expression (12) and taking the temporal Fourier 
transform of equation (17) at the frequency lo, one obtains 

ap(cu) = (2'7r)‘* ^ [a< 3 (w) + ^^(t^)] x 

<1 

E /dJ ■ (18) 

m 

We define the response matrix of the system M as 

Mp,(c.) = (27r)‘'E/ 

m 

where one must note that the response matrix depends 
only on the initial equilibrium state of the disc, since 
dFo/dJ evolves only on the secular scale, the perturbing 
and self-gravitating potentials are absent, and the basis el¬ 
ements are chosen once for all. Finally, in order to 

shorten the notations, the amplitudes of the self and exte¬ 
rior potentials are defined as a(t) = (ai(t),..., ap(t),...) and 
b{t) = (6i(t),..., bp{t), ...). Thanks to these notations, one can 
simplify equation (18), and rewrite it under the form 

a(cj) + b(a;) = [l-M(a;)]~’‘-b(a;). (20) 

One should note that the matrix [I —M] is invertible only if 
the self-gravitating system is linearly stable so that all the 
eigenvalues of M are assumed to be strictly smaller than 1 
for all values of uj. 


2.5 Diffusion coefficients 

The amplification relation (20) corresponds to the short 
timescale (dynamical) response of the system, driven by the 
evolution equation (7). We will now describe the impact of 
these solutions on the long timescale diffusion equation given 
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by equation (6). Starting from equation (6), one has to eval¬ 
uate an expression of the form 


(27r) 






( 21 ) 


(27r) 




Here, only terms with mi = —m 2 are different from 0. Using 
equation (10) and the fact that ip-m = '4>m, we can finally 
rewrite the diffusion equation (6) under the form 


dFo 

dt 



m 


D, 


, n dFo' 


( 22 ) 


where the anisotropic diffusion coefficients Dm{J) are given 
by 

-Dm(J,t) = [V’^*(J,t)+l/’^*(J,t)] X 

/ t ( 23 ) 

[V>)).(J,r)+^;,(J,r)] . 

-00 


Note that equation (22) can be re-arranged as 


dFo 

dt 


d 

dj' 




, with D(J) m (J) m(g)m, 


an anisotropic tensor diffusion matrix. Using the basis de¬ 
composition introduced in equation (15), the diffusion coef¬ 
ficients from equation (23) take the form 

P,(l 

f dr [ap(r)-|-6p(r)] . (24) 

J —00 


Expressing the temporal coefficients ap{t) and hp{t) via their 
Fourier transforms, we obtain 


. iut . , 

e X 


^dcj'[ap-|-bp] (cu') . (25) 


dr e 


— im-d(t—T) j 


The amplification relation (20) allows us to rewrite equa¬ 
tion (25) as 


^ PyQ Pi ,91 

|dc.e*“‘[[l-M(^)];;J‘fe;(a,)x 


(26) 


/ dre 

' —00 


— im-d(t—T) i 


dPe [I-M(cu')] I bp^ {P) ■ 


2.6 Statistical expectation 

The final stage of the derivation is to introduce the statistics 
of the external perturbations. Indeed, our previous calcula¬ 
tion corresponds to the response of the system to a given 
particular perturbation history: 1 1 —> b{t). Let us now de¬ 
note the ensemble average operation on such different real¬ 
izations as (.). As the global underlying background gravi¬ 
tational potential is assumed to be stationary, the mapping 
(x,v) f-,- (0,J) remains the same for the different realiza¬ 
tions, so that the operations of derivation or integration with 


respect to 9 and J commute with the ensemble average. The 
diffusion equation (22), when ensemble averaged, takes the 
form 


d{Fo) 

dt 



m 




dFoV 

' dJ /. 


(27) 


A priori, the gradient dFo/dJ, cannot be taken out of the 
ensemble average operation. However, we intend to describe 
the effect of an averaged fluctuation on a given Fo repre¬ 
senting a mean disc. Then, one may assume that the quasi¬ 
stationary distribution function To, its gradients and there¬ 
fore the response matrix M do not change significantly from 
one realization to another, so that they can be taken out of 
the ensemble average operation. This means that we assume 
there exists a mean response for the secular distribution. 
To = (To), de-correlated from the perturbations, so that we 
have {Dm{J) m-dFo/dJ) = {Dm{J)) m-dFo/dJ . We also 
suppose that the time evolution of the exterior perturbing 
potential is stationary and therefore introduce the corre¬ 
sponding temporal autocorrelation function defined as 

Cki{ti—t2) = {bk{ti) b*{t2)) , (28) 


where the exterior perturbation is also assumed to be of 
zero mean. The autocorrelation C connects the temporal 
coefhcients b, whereas the diffusion coefficients from equa¬ 
tion (26) involve the Fourier transformed ones b, so that one 
needs to compute (&fe(w) 6*(w4)- One can straightforwardly 
show that 

(bk{uj)bi{uj')) = 2tv 5o{u!-uj')Cki{oj) , (29) 


where C is the temporal Fourier transform of the autocor¬ 
relation of the external potential. Using this result in the 
ensemble averaged expression (26) yields 


{Dm{J,t)) = f do; f dr'( 

^ ^ J J — 00 


—i (cu — m • £1) T ^ 


[I-M]"bc-[I-M]"^l (tu), (30) 

J PQ 


where we performed the change of variables t' = r — t. 
One should note that when ensemble averaged, the diffu¬ 
sion coefficients are (explicitly) independent of P. In or¬ 
der to shorten temporarily the notations, we introduce 
L = [I— M]“^-C- [I— Mj”^. In equation (30), one has to 
evaluate an expression of the form 


1 


da; L(a;) / dr^e 


27r 


‘-X 


+ 00 

dai 


L(^) 

Ld — m fl ’ 


(31) 


where in the integration over t' we only kept the term for 
r' = 0, by adding an imaginary part to the frequency 02 so 
that uj = u!+i0'^, ensuring the convergence for t'^—oo. 
The remaining integral over oj can be evaluated using 
Plemelj formula 


1 

X ± i0+ 



=F iTTSo{x) , 


(32) 


^ though they depend on the secular timescale via the variation 
of F in M. 
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where V denotes Cauchy principal value. Therefore, equa¬ 
tion (31) becomes 

(31)(x^p/dX +lL(m.n). (33) 

Ztt Juj — m ll Z 


The final step of the derivation is to show that the prin¬ 
cipal value term present in equation (33) has no impact 
on the secular diffusion. Indeed, using the expression (23) 
of the diffusion coefficients, one can show that they satisfy 
D-rn{J) = D’^{J). As a consequence, as we are summing 
on all the modes m, the diffusion equation (22) may be 
rewritten under the form 


dFo 

dt 


E’ 


dj 


Re [Dm,{J)] m 


dFo 

dJ 


(34) 


From equations (19) and (28), we know that the response 
matrix and the autocorrelation matrix are hermitian so that 
M*=M* and C*=C*. As a consequence, the matrix L 
is also hermitian. Since Re(D m ) = {D m + 49m)/2, starting 
from equation (33), we immediately obtain 

(Re[Z?^(J)]) 

P,(l 

C(m-fi), (35) 

L J pq 

SO that the full secular diffusion equation takes the form 


dFo 

dt 


m *- p^q 

C-[I-M]"^l {m-n) . 

L Jpg _ 


(36) 


Equation (36) is the main result of this section. In Ap¬ 
pendix A, we present an alternative derivation of these diffu¬ 
sion coefficients based on Hamilton’s equations from which 
we recover the exact same diffusion equation. The 1/2 fac¬ 
tor recovered via these two complementary approaches was 
skipped in the calculation presented in Pichon & Aubert 
(2006), because of an error in the bounds of half-temporal in¬ 
tegrations, similar to the one present in equation (31). This 
derivation is valid in any dimensions, provided the under¬ 
lying system is integrable. One may also note that in the 
homogeneous limit, equation (36) reduces to the secular dif¬ 
fusion equation obtained in Chavanis (2012b, 2013). In the 
next section we will restrict ourselves to 2D conhgurations 
and make further assumptions in order to simplify equa¬ 
tion (35) into a one dimensional quadrature. 


3 THIN TEPID DISCS AND WKB LIMIT 

One difficulty for the implementation of the secular dif¬ 
fusion equation (36) is to simultaneously have an explicit 
mapping {x,v)i-^{0,J) to the angle-actions coordinates, 
and be able to evaluate the diffusion coefficients given by 
equation (35), which require to invert the response matrix 
[I —M]. In order to deal with the non-locality of Poisson’s 
equation, we also have to explicitly introduce potential basis 
elements, as in equation (14), to compute the response 
matrix from equation (19). To ease these calculations in a 
2D axisymmetric disc, one may rely on the WKB assump¬ 
tion (Liouville 1837; Toomre 1964; Kalnajs 1965; Lin & Shu 


1966; Palmer et al. 1989), which assumes that the perturba¬ 
tions and self-responses will take the form of tightly wound 
spirals, which in turn allows us to write Poisson’s equation 
locally. Considering only such perturbations sums up to in¬ 
troducing basis elements with specihc properties as detailed 
later on. 


3.1 Epicyclic approximation and isothermal DE 


In order to explicitly build up a mapping (x, v) !->■ (6, J) for 
an axisymmetric disc, we assume that the disc is sufficiently 
cold and therefore rely on the so-called epicyclic approxima¬ 
tion. 

The natural coordinates for an axisymmetric galactic 
disc are the polar coordinates {R,4>), with their associated 
momenta {pR,p^). Within such coordinates, the stationary 
Hamiltonian of the system reads 


Ho{R,<t>,PR,P4’) 


1 

2 




+ V’o(R), 


(37) 


where V’o is the axisymmetric stationary background poten¬ 
tial within the disc. The hrst action of the system is the 
angular momentum dehned as 


Ja} — Lz — 


1 


A(f)PA, =Pa> = R (j>. 


(38) 


As soon as the value of J^, is imposed, one obtains a new 
equation of motion for the R variable given by 

d'lpeB 


R = -- 


dR 


where the effective potential is defined as 


■tpesiR) = tpo{R) + 


04 , 

2^ ■ 


(39) 


(40) 


The main idea behind the epicyclic approximation is to ap¬ 
proximate the radial motion as an harmonic oscillation. For 
a given value of we dehne implicitly the guiding radius 
Rg as 


9V’eff 

dR 


Rg 


d^po JI 


(41) 


so that Rg{J 4 ,) corresponds to the radius for which stars 
with an angular momentum equal to J 4 , evolve on circular 
orbits. For a stationary potential, the mapping between Rg 
and J 4 , is bijective and unambiguous (up to the sign of Jf). 
We define the azimuthal frequency Q.{Rg) and the epicyclic 
frequency K{Rg) as 


n^Rg) 

K\Rg) 


1 dlpo 

'Wg~^ 

9^V’eff 
dR^ , 


t2 

'Jch 


Rj 


d^tpo 


-u 

R'^’ 

R„ “9 


(42) 


A Taylor expansion at first order near Rg of equa¬ 
tion (39) shows that R satisfies the differential equation 
R = —K^{R—Rg), which is the evolution equation of an har¬ 
monic oscillator centered on Rg. We introduce the amplitude 
A of the radial oscillations and define the radial action Jr 
as 




1 

2^ 


dRpR = 2 '^^ 


2 


(43) 
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The case Jr = 0 corresponds to circular orbits. The larger 
Jr, the wider the radial oscillations of the star. One should 
note that within the epicyclic approximation, the two in¬ 
trinsic frequencies and n are only function of the angular 
momentum J^ and do not depend on the radial action Jr- 
Finally, one can show (Lynden-Bell & Kalnajs 1972; Palmer 
1994; Binney & Tremaine 2008) that the mapping between 
{R, 4>,Pr,P4>) and {Or, 6 ^, Jr, J 4 ,) takes at first order the form 


R = Rg+Acos{9R), 

A, a 211 . 

0 = 9 ^ -— sm(6»H). 

K Rg 


(44) 


Within this approximation, one can easily parametrize plau¬ 
sible stationary distribution functions for a galactic disc, 
defined as functions of the actions (J,/,, Jr). Indeed, we sup¬ 
pose that the stationary distribution Fq of the disc is a 
Schwarzschild distribution function (or locally isothermal) 
given by 


Fo{Rg, Jr) 


n{Rg)EiRg) 

TT R{Rg) (T^{Rg) 


^{Rg) Jr 
a^{Rg) 


(45) 


where Ti{Rg) is the surface density of the disc and a^{Rg), 
which depends on the position in the disc, encodes the typ¬ 
ical radial velocity dispersion of the stars at a given radius. 
The larger the hotter the disc and the more stable it is. 


3.2 The WKB basis 

In order to use a WKB approach in the secular diffusion 
equation (22), one needs to introduce explicitly a basis of 
density-potentials from which the WKB hypothesis will fol¬ 
low. 


3.2.1 Definition of the basis elements 

We define on the plane the following basis of potential func¬ 
tions, well suited to represent tightly wound spirals 

^[H,kr,Ro\ 4))=A ^ ( 45 ) 



Figure 1. Two WKB basis elements. Each Gaussian is centered 
around a radius Rq . The typical extension of the Gaussian is given 
by the decoupling scale cr, and they are modulated at the radial 
frequency kr. 


3.2.2 Assoeiated surface density elements 

Equation (14) requires to have a biorthogonal potential ba¬ 
sis. We now determine the surface density basis elements as¬ 
sociated to the potentials from equation (46). For 2D razor 
thin discs, we are in the presence of a discontinuity leading 
to a surface density S(i?, 0). In order to satisfy Poisson’s 
Equation, we extend our potential to the g-axis via the fol¬ 
lowing Ansatz 


Injecting such an expression into Poisson’s Equation in vac¬ 
uum = 0, we obtain after some algebra 


R — Ro 1 R — Rq 




kr R 

(ak, 


1 , kl 

R — Rq 1 

2 - 



{akr)^ (krR)^ 

Cr2 kr 




+ 

(49) 


In order to obtain a simple expression for the surface density 
basis elements, we introduce additional WKB-like assump¬ 
tions, so that the terms appearing in equation (49) are all 
negligible in front of 1. First of all, we assume that the spirals 
are tightly wound so that we have 


where the functions Br^ (R) are of the form 


krR ;?> 1. 


(50) 


BRoiR) 




(R-Ro)^ 

2cr2 


(47) 


To our knowledge, this is the first time that such tightly 
wound basis elements have been introduced in the context 
of discs secular dynamics. The central radius Rq is the ra¬ 
dius around which the Gaussian Brq is centered, is an az¬ 
imuthal number representing the angular component of the 
basis elements, and kr corresponds to the radial frequency 
of the potential element. Here cr is a scale-separation pa¬ 
rameter ensuring the biorthogonality of the basis elements 
(see below). The radial dependance of the basis elements is 
illustrated in figure 1. The amplitude A is not yet defined 
and will be chosen in order to guarantee the correct normal¬ 
ization of the basis. The unusual normalization of Br^ will 
ensure that A is independent of cr and will also allow us to 
naturally introduce Dirac deltas 5r,{R—Ro) in some of the 
next calculations. 


Moreover, introducing the typical size Flays of the system, 
we add the supplementary constraint 

kr<J^ —— . (51) 

cr 

In this limit, assuming that k^ remains of the order of the 
unity, equation (49) becomes 



Therefore, we conclude that within the WKB approxima¬ 
tion, the extended 3D potential can be written as 

'tP^'°<i>’'°-’^«\R,(j,,z) ='tP^'°*’'°-’^«]{R,<f)e-^’-\^\, (53) 

where we added absolute value on « in order to respect the 
boundaries conditions of the potential at 2 = ±cx), where the 
potential has to tend to 0. Such a potential introduces a dis¬ 
continuity of d'4>/dz at the plane 2 = 0, consistent with the 
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given surface density. Gauss theorem for the discontinuities 
at a plane may be written as 


E(7?,0) 


1 


y dS 

r^0+ OZ 


lim 
2—^0~ 


dij) 

dz 


(54) 


We immediately conclude that the surface density associated 
to a given potential element jg given by 




i[k^,kr,Ro] 

2tvG ^ 




(55) 


With these assumptions, we can ensure that we must neces¬ 
sarily have — kt- and 7 ?q = Rq, in order to have 

a non-zero term. The last step is to explicitly calculate the 
amplitude A of the basis elements. Indeed, starting from 
equation (56), we have the condition 


— ^ [dR R exp 

G V^cr j 


{R-Rof 

n-2 


(61) 


which may be rewritten as 



1 — ^0 

l-ferf 

Rq 

1 a _f? 2 / 2 

+ e «o/- 

3.2.3 Biorthogonality condition 

~ G 2 

- 

(7 

■/tt Ro 


(62) 


The next step of the definition of the WKB basis is to ensure 
that this basis is biorthogonal as in equation (14). Indeed, 
it has to satisfy the property 


fcl Ul 


dRRd(j) . 


The r.h.s. of this expression becomes 


Using the assumptions made in equation (60), we imme¬ 
diately conclude that 7?o/u ^ 1, so that erf [i?o/o'] ~ 1 and 
(j/{^/ tvRq) exp[—i?o/cr^] <C 1. We therefore finally obtain the 
expression of the amplitude of the basis potentials as 




/ 


G 

\kr \ Ro 


(63) 


jfl 

2tvG 


Ap An 





X 


(56) 




r (R-Kf] 


r (R-Rif] 

2cr2 

exp 

2cr2 


The integration on (j> is straightforward and gives a term 
fcP 

equal to 2tv5 A . In order to be able to integrate on R, we 

now need to impose new WKB-like assumptions to justify 
the biorthogonality of the basis. We introduce the spatial 
Fourier transform R with respect to R, and the difference 
of the two radial frequencies Akr = k^ — k^. Equation (56) 
requires to integrate an expression of the form 


= I dk' R[By,pf{k')R[By,.f{Akr-k') 
o^jdk'ej 


(57) 


2/a 


2/a 


-i(Rlk'+R%[/\kr-k']) 


3 . 2.4 Fourier development in angles 

The diffusion equation involves terms of the form ipl/f (J), 
that we will now evaluate for the WKB basis, equation (46). 
Using the epicyclic mapping from equation (44), we need to 
estimate 

fde^ x (64) 

G-r) J J 

The integration on 0,^ is straightforward and is equal to 
2-17Sf. Looking at the dependence in On within the com¬ 
plex exponential, we can write 

krAcos{9R)-k^ — ^ sm{9R) = H^J/kr) fiin{9R+9f) , 

where we have defined 


where we used the property that the Fourier transform of a 
product is given by the convolution of the Fourier transforms 
and that the Fourier transform of a Gaussian of spread a is 
a Gaussian of spread 1/a. In expression (57), we note that if 
\k'\ ^ 1/ct or \Akr — k'\ ^ 1/a, then the product of the two 
terms can be considered to be negligible. We will therefore 
suppose that one of the conditions 

Akr S> — or Akr = 0 , (58) 

a 

holds. Under this assumption, the term is non-zero only for 
k/ = k/. Finally, it remains to prove that non zero terms are 
only obtained when Rq = Rf The peaks of the two Gaus- 
sians in equation (56) can be considered as sharp and sepa¬ 
rated if Ai?o = Rq—Rq satisfies the condition 

Ai?o ^ (T or ARo — 0. (59) 

To sum up our assumptions so far, we should consider peak- 
radius Ro, spread a and radial frequencies kr such that 

ARo > cr > . (60) 

To these conditions, one must also add the constraints ob¬ 
tained in equations (50) and (51) via Poisson’s equation. 


Hkfkr)=Ajk/ + kl 


20 . 


txRn 


9% = tan ^ 


K krRg 

20 k^ 

(65) 


Thanks to our WKB assumptions, an approximation of the 
amplitude term Hkfkr) and the phase-shift 9o is possible. 
Indeed, both of these terms involve an expression of the 
form 2k^xO/RXl/(krRg). Yet, we made the assumption 
that krRg 1. Moreover, we know that for typical galaxies 
(Binney & Tremaine 2008). Assuming that 
k^ is of the order of unity, we obtain the approximations 


Hkfkr)-A\kr\-^l^\kr\ ', " J • (66) 

We also supposed that the radial oscillations are small, so 
that the epicyclic amplitude satisfies A Rg. Thanks to 
this assumption, we may replace BRg{Rg +Acos{9Ry) by 
RnoiRg), keeping the dependence on A only in the complex 
exponential. This is a crucial step to be able to integrate 
explicitly on 9r. We also introduce the Bessel functions of 
the first kind which satisfy the property 

^iH^^ikr) sin(flK + 0?;) [Hkfkr)] + . (67) 
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It is then possible to perform explicitly the integration on Or 
in equation (64), which is equal to 2n 5™'", so that only one 
Bessel function remains. We finally obtain the expression of 
the Fourier transform in angles of the basis elements 


3.3 Estimation of the response matrix 


Using the explicit WKB potential basis introduced in equa¬ 
tion (46), one can now estimate the matrix response from ex¬ 
pression (19). The approximation obtained in equation (66) 
allows us to simplify the phase-shift terms, so that equa¬ 
tion (19) becomes 








(jj — mfl 


{H^AK)) Brp (r,) Br, {r,) . 


One should note that this expression is similar to equa¬ 
tion (56). Indeed, using the assumptions from equation (60), 
we are able to ensure that only the diagonal coefficients of 
the response matrix are different from 0. First of all, the 
azimuthal Kronecker symbols impose that There 

is however a slight complication in the calculation because 
of the presence of additional terms depending on Rg. In or¬ 
der to sketch the proof of this statement, we introduce the 
function 


h(Rg)^ 



m-dFo/dJ 

dRg 

uj — m fl ' 


'ApAqAm 




Jm 


kl 


(70) 


This function captures all the additional Rg dependence 
appearing in equation (69). Using the change of variables 
1-^ Rg, the integral on which has to be evaluated in 
equation (69) takes the form 


IdRgh{Rg ) exp 


2a2 


2cr2 


(71) 


Using the assumption from equation (59) relative to the pos¬ 
sible values of A7?o in the WKB basis, one can note that the 
product of the two Gaussians in Rg imposes Rq = Rq in or¬ 
der to have a non-zero contribution. The previous expression 
then becomes 


(71) oc JdRg h{Rg) exp 


{Rs-Rof 


(72) 


Using the same argument as in equation (57), we can rewrite 
the Fourier transform in R as the convolution of two radial 
Fourier transforms so that it becomes 


(71) cx 


dk' R[h]{k') exp 


{Akr-k'f 

^ _ 


(73) 


We now use equation (58) relative to the possible values of 
Akr in the WKB basis. If we suppose that Akr A 0) the 
width of the Gaussian from equation (73) imposes that the 
integration will only probe the contribution of J-[h] in the 
neigborhood of k' ~ Ak^ ':$> 1/a. We assume that the radial 
Fourier transform of the function h is such that it is mainly 
focused in the frequency region |fc| < l/a, meaning that for 
a typical galactic disc, the main frequencies of radial varia¬ 
tions of h are inferior to l/a. Under this assumption of slow 


radial variation within the disc, one can see that non-zero 
contributions can only be obtained for Akr = kr — kr = 0. 
As a consequence, we have shown that within the WKB 
approximation the response matrix is diagonal. For these 
diagonal coefficients, it only remains to evaluate explicitly 
the integrals over J,j, and Jr in order to obtain the expres¬ 
sion of the response matrix eigenvalues. This calculation is 
presented in Appendix B. Within the assumption that the 
galactic disc is tepid, the eigenvalues of the response matrix 
take the form 


M 


[kl,kl,Ro\,[kl,k'i,Ro\ 


{uj)=5 


k^r^l 27rGS \ kr\ 


R{s,x), (74) 


where x s-nd s are respectively defined in equations (B6) 
and (BIO) and F{s, x) is the reduction factor introduced 
in equation (Bll). This amplification eigenvalue is in full 
agreement with the seminal works from Kalnajs (1965) 
and Lin & Shu (1966), which independently derived the 
WKB dispersion relation for stellar discs.® 


3.4 Estimation of the diffusion coefficients 

The expression (35) of the diffusion coefficients shows that 
the diffusion coefficients require the evaluation of [I—M]“^. 
In order to simplify the notations, we will denote our poten¬ 
tial basis with only one index so that 

Equation (74) shows that the response matrix is diagonal in 
the WKB approximation. We therefore introduce the eigen¬ 
values of M as 

Ap = Mpp . (76) 

The matrix [I—M]“^ is then diagonal and reads 

[I-M]-i = 5® ^ . (77) 

Thanks to these diagonal coefficients, the expression of the 
diffusion coefficients from equation (35) becomes 

= ^-*Y^Y^Cp,(m.fl). (78) 

P.9 ^ 

At this stage, we use the property from equation (29) to 
rewrite Cpg as a function of the basis coefficients bp and 
bg. Remembering that the basis elements ifm and the ma¬ 
trix eigenvalues Ap do not change from one realization to 
another, one can rewrite equation (78) under the form 

Drr.iJ)=(^ 

j^j^bp{m-n)b*g{uj')j. (79) 

It is important here to note that the eigenvalues Ap, Xq and 
the basis coefficient bp are both evaluated at the intrinsic 


® For nice introductions to the WKB dispersion relation in stellar 
discs, see section 6.2.2 of Binney & Tremaine (2008) and section 
1.4.2 of Binney (2013). 
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frequency m Cl, whereas b* has to be evaluated at the in¬ 
tegrated frequency lo'. In the upcoming calculations, in or¬ 
der to shorten the notations, when obvious, the frequencies 
of evaluation will not be written. Using the expressions of 
the basis elements in the WKB approximation from equa¬ 
tion (68), we can write 


‘ if; Y. \ 


G 




Jm 




Jm 




JRoikt-kl) I 1 


: exp 


r [Rs-Rif] 


r [Ra-R^of] 

2cr2 

exp 

2(t2 


bpb; 


(80) 


where we already got rid of the sum over and since 
the Fourier transform of the WKB basis elements from equa¬ 
tion (68) imposes to have 


m^ = kl = kl. (81) 

In the expression (80), we also neglected the phase terms in 
nirOff and simplified the value at which the Bessel functions 
have to be evaluated using the approximation introduced in 
equation (66). 

In order to obtain an expression independent from the 
choice of the basis (i.e. the precise value of a), we will now 
replace the coefficients bp by their expressions in terms of 
the true exterior potential function i/)®, which is completely 
independent of the choice of the basis. As the potential basis 
in the WKB approximation is bi-orthogonal, the temporal 
Fourier transform of the basis coefficients is given by 

bp{oj) = -Jd^X ip‘^(x,Lo), (82) 

where the hat ^ corresponds to the temporal Fourier trans¬ 
form defined in equation (11). Using the expression of the 
surface density basis from equation (55), we obtain 

bp=ldRRld^^A {R) r (R, 4>) ■ 


The integration on (j> is straightforward and leads to a term 
equal to 27r '!/’®(.p (R), where the presence of the index k^ cor¬ 
responds to the Fourier transform with respect to the phys¬ 
ical angle (j>, using the same conventions as in equation (8). 
We may now write 



dRR exp 




2(t2 


e-^“?^|p(7?). (83) 


This integration should be interpreted as the radial Fourier 
transform at the frequency fcjl of the exterior potential in the 
region close to Rq. Since the integrand contains a Gaussian 
in R of spread a, we may take the term in R out of the 
integral and consider it to be equal to Rq. We now define 
the local Fourier transform of the exterior potential on a 
restricted region of radius (Gabor 1946) as 


i’%^,kr[Ro] 


1 

2711 


dRiplJR) exp 


[R-Rof 

2cr2 




This definition is motivated by the fact that if we consider 
the case of an uniform perturbation = 1, then its local 
Fourier transform is independent of Rq. One may note that 
this definition is not independent of the decoupling scale a, 


but as we will see later on, it is the relevant quantity in order 
to obtain diffusion coefficients independent of this ad hoc 
parameter. Thanks to this definition, the basis coefficients 
from equation (83) become 



271 

(tTCT^)!/! 



(85) 


We recall the notation used for the exterior potential in the 
previous expression. The index k^ corresponds to the az¬ 
imuthal Fourier transform with respect to the physical angle 
(j), and the index fcj? corresponds to the local radial Fourier 
transform with respect to the physical radius R in the neigh¬ 
borhood of Rq, as defined in equation (84). The diffusion 
coefficients from equation (80) are then given by 


DMJ)=(±fdc7'J2jr 


J„P 


J™ 




i(Rg-RP)kP i(R„-R7)k7 1 1 


1 Ap 1 \q 


27V 

^exp 


r [Rg-R^Qf] 


r (Ra-Rlf] 

2cr2 

exp 

2a2 


'4’‘^ms,k?[Ro]S'^m^,kfRo] ) ' 


( 86 ) 


One can note that the gravitational constant G has disap¬ 
peared, since the dependence on the strength of the gravity 
is now hidden in the units of i)’®. One should also recall 
that the previous expression has to be evaluated at the res¬ 
onance frequency, so that lo = m Cl, except for k?.[Ro] 
which has to be evaluated at the frequency ui'. The main 
step of the simplification is now to replace the discrete sums 
on the basis index kr, kt, Rq and Rq by continuous inte¬ 
grals. One should indeed now recall that our potential basis 
elements are made of three different index. Here k^ is a dis¬ 
crete index which must necessarily be equal to m,^, so that 
it is absent from the sums, kr is a continuous index, whose 
value has to belong to ]1 /ct ; ...[, because of the approxima¬ 
tions made in equation (60), and finally Rq whose values 
belong to ]a; ...]. We must also comply with the two as¬ 
sumptions (58) and (59) about the distance Akr and ARq 
between two consecutive elements of the basis. In order to 
get rid of the sum over the discrete index, we will use Rie- 
mann formula f{x)Axc:^ fdxf(x), with Aa; controlling 
the distance between two consecutive elements. The depen¬ 
dences with the two radial frequencies kr, kr and the two 
radii Rq, Rq are such that the sums on the index p can be 
completely disentangled from the sums on the index q. In 
order to emphasize the gist of the calculation, the diffusion 
coefficients from equation (86) may be written under the 
form 

DrrviJ) = (^-^jdLj'g{m-n)g*{uj'^ , (87) 

where g{ijj) is defined as 

g{uj) = Q{Rg-Rl). (88) 

k^,Rl 

In equation (88), gs{kr, Rq,lo) encompasses all the slow de¬ 
pendences of the diffusion coefficients with respect to the 
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position Ro and the radial frequency kr so that 

1 


Qsik^^ Rq: — f/rn 




1-A. 


(89) 


and Q{Rg — Rq) is a normalized Ganssian given by 


g{Rg-Rl) = 




exp 


{Rg-Rlf 


2a2 


(90) 


In the discrete sum from equation (88), the basis elements 
are separated by constant step distances A_Ro and Afcr. We 
suppose that generally fc)? and Rq are given by 


= UkAkr , 

Rq = Rg + TItARq . 


(91) 


where Uk is a strictly positive integer and Ur is an integer 
that can be both positive and negative. One can note in 
equation (88) the presence of a rapidly evolving complex 
exponential which may cancel out the diffusion coefficients 
if the basis step distances are not chosen carefully. Injecting 
the dependences from equation (91) in the complex expo¬ 
nential from equation (88), one can see that we have to sum 
terms of the form 

exp(i{Rg — {Rg+nrARo))nkAk^ = exp(—inrnkARoAk^ . 

As a consequence, since TirTik is an integer, in order to have 
no contributions from the complex exponential term, one 
has to choose step distances so that 


Ai?o Afcr = Stt . 


(92) 


This choice of step distances, imposed by the complex expo¬ 
nential term, corresponds to a critical sampling (Daubechies 
1990), which allows us when performing the change to con¬ 
tinuous expression in equation (88) to leave out the complex 
exponential. This transformation is a subtle stage of the cal¬ 
culation, since we require our step distances ARo and Afcr 
to be simultaneously large to comply with the WKB con¬ 
straints from equation (60) and small to allow the use of 
Riemann sum formula. As the radial Gaussian g{Rg — Rg) 
is sufficiently peaked and correctly normalized, one can re¬ 
place it by So{Rg — Ro)- The integration on Rq can then be 
immediately performed to obtain 

g(uj) = J dfc(! fls(fcr, Rg, u) . (93) 


Using this result in equation (87), we hnally obtain the ex¬ 
pression of the diffusion coefficients as 


Dm(j) 



dio'Jdk^Jrr, 



1 

l--^fc? 



[Rg] ^ 



1 




In this expression, all the radial functions have to be evalu¬ 
ated at the position Rg and at the temporal frequency m Cl, 
except for which is evaluated at the frequency 

oj'. The eigenvalues Afcr are given by equation (74) and read 


Afcr [Rg,m-fl] 


27rGE|fcr| 

k2(1-s2) 


ns,x)- 


(95) 


Note that, as requested, in equation (94), all the dependen¬ 
cies in a have disappeared, so that the value of these diffu¬ 
sion coefficients is independent of the precise choice of the 


WKB basis. One can finally introduce the autocorrelation 
of the external pertubation C,/, as 


C,jj[m,l,,UJ,kr,kr, Rg] = 


(96) 



,kP[Rg,^]'<Pn 



so that the expression (94) of the diffusion coefficients takes 
the form 


Dm{J) = Jdk^rJmr 

[PP kP 

y K R 

jdKJmr 

,fppjp 

y K R 


l“-^fc; 

1 


1“ Afc 


- X (97) 

■G,/,[m,/,,m-I2,fcr,fcr,Rg]. 


We may hnally asume that the external perturbations are 
spatially quasi-stationary, so that we have 


[-Ri, ti] ]R2S2])= (98) 

C[m,p,ti—t 2 , Ri —R 2 , {Ri+R2)/2] , 


where the dependence of C with respect to {Ri+R2)/‘2 is 
supposed to be weak. As demonstrated in Appendix G, one 
can then show that 

(^m,^,fcl[^9.^l] V''C^,fc2[R9-^2]) = (99) 

27r So{aJi-aJ2) So{kl.—kr)C kl.,Rg]. 


Using this autocorrelation function diagonalized both in cu 
and fcr, the expression of the diffusion coefficients from equa¬ 
tion (94) Hnally takes the form 


7 

O rr\ 


dkr 


[l-Afcr]' 


C[m^,m-fl,kr,Rg]. (100) 


Equation (100) is the main result of this section. The corre¬ 
sponding anisotropic tensor diffusion coefficient reads 


D = 




/ dkr 



[l-Afc.]" 


C[mrl,,Tn-fl,kr,Rg] . 


One may sometimes simplify further equation (100) when 
the function kr 1 —>• Afc,, is a sharp function reaching a maxi¬ 
mum value Amax(Rg,l^ = mfl), for kr = kraBxiRg , aj), with 
a characteristic spread given by Akx. Under this assumption 
of so-called small denominators, the previous expression of 
the diffusion coefficients can be approximated as 



Dm{J) = Akx—^ - 

[1 '^fcmax 


H [nifj), Tjrt • 12, fcmax, Rg] • (101) 


One should note here that the autocorrelation of the exter¬ 
nal perturbation C, which sources the diffusion coefficients 
Dm (J) depends on four different parameters: the azimuthal 
wave number the location in the disc via Rg, the radial 
frequency fcmax of the most amplihed tightly-wound spiral at 
this position, and Hnally the local intrinsic frequency m fl. 


4 DISCUSSION AND CONCLUSION 

Starting from Boltzmann’s collisionless equation expressed 
in angle-actions coordinates and relying on a timescale de¬ 
coupling, we derived in equation (36) a diffusion equa- 
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tion describing the long-term evolution of a perturbed self- 
gravitating collisionless system^. This general formalism is 
appropriate to capture the nature of a collisionless system 
(via its natural frequencies and susceptibility) as well as its 
nurture via the structure of the power-spectrum of the ex¬ 
ternal perturbations. Hence, it yields the ideal framework in 
which to study the long-term evolution of such system. 

When applying this Fokker-Planck diffusion equation 
to an infinitely thin galactic disc, we used two main ap¬ 
proximations. We first assumed the disc to be tepid. Having 
orbits with small radial oscillations justified the use of the 
epicyclic approximation, allowing us to explicitly build up in 
equation (44) a mapping between the physical coordinates 
{x,v) and the angle-actions coordinates (0, J). Another im¬ 
portant consequence of the epicyclic development is to allow 
for a direct determination of the local frequencies of the sys¬ 
tem Q and k, as in equation (42). Being able to localize the 
resonances is crucial in this formalism, since the diffusion 
coefficients from equation (36) show that both the suscepti¬ 
bility of the system via [I —M] and the external perturbing 
power spectrum via C have to be evaluated at the intrinsic 
frequency m fl. The second approximation involves an ex¬ 
plicit WKB basis introduced in equation (46). It allowed us 
to obtain in equation (74) a diagonal response matrix, as if 
gravity was only local. Thanks to the assumption of radial 
decoupling, the WKB approximation led to equation (100), a 
simple quadrature for the diffusion coefficients, with which it 
is straightforward to identify the physically relevant modes. 
Such simplification provides useful insight into the physical 
processes at work, e.g. the relevant resonances, their loci and 
their relative strengths. 

The formalism of secular resonant dressed orbital dif¬ 
fusion and its WKB limit is implemented in the companion 
paper (Fouvry & Pichon 2014, paper H, submitted) to re¬ 
cover the formation of resonant ridges in action-space when 
an isolated stellar Mestel disc (Mestel 1963) is left evolving 
for hundreds of dynamical times. The development of such 
ridges has been shown to originate from a resonant mono¬ 
dimensional diffusion, specifically enhanced in restricted lo¬ 
cations in the disc. It captures the respective roles and im¬ 
portances of various parameters of the system. Indeed, pa¬ 
per H illustrates on an example that the self-gravity of the 
disc (via the amplification eigenvalues A), its susceptibility 
(via the anistropic diffusion coefficients Dm{J)), its inho¬ 
mogeneity (via the gradients dFo/dJ), its temperature (via 
cr^), its physical structure (via the introduction of taper¬ 
ing functions representing resp. the bulge and the outer 
edge of the disc), and the detail of the source of pertur¬ 
bations (via the power spectrum of i/>®), all contributes non- 
negligibly to the appearance of resonant ridges. Such fea¬ 
tures have been observed both in numerical experiments 
(Sellwood 2012) and in the Solar neighborhood (Wielen 
1977; Dehnen 1998; Nordstrom et al. 2004; Famaey et al. 
2005; Aumer & Binney 2009; McMillan 2011). 

The WKB assumption can also be used to study the col- 
lisional evolution of a self-gravitating disc containing a finite 
number of substructures. Indeed, in Fouvry et al. (2014a, in 
prep.), the same local WKB approach will be applied to the 

^ Appendix A also shows how this diffusion equation is obtained 
via a different route involving Hamilton’s equations. 


Lenard-Balescu non-linear equation (Balescu 1960; Lenard 
1960; Weinberg 1998; Heyvaerts 2010; Chavanis 2012a), 
which accounts for self-driven orbital secular diffusion of a 
self-gravitating system induced by an intrinsic shot noise 
due to the discreteness of the system. Possible cases of ap¬ 
plications of this approach are the secular diffusion of giant 
molecular clouds in galactic disc, the secular migration of 
planetesimals in proto-planetary discs, or even the long-term 
evolution of population of stars within the Galactic center. 

For self-gravitating systems which do not take the form 
of an infinitely thin disc, for which the epicyclic approx¬ 
imation and the WKB assumption may be relevant, the 
formalism of secular forcing can still be used. The diffu¬ 
sion equation (36) could for instance describe the secu¬ 
lar diffusion of dark matter cusps in galactic centers in¬ 
duced by perturbations from stochastic feedback processes 
originating from the baryonic disc (Fouvry et al. 2014b, 
in prep.). Given a detailed characterization of the per¬ 
turbations induced by e.g. the cosmic environment, one 
could also study their long-term effects on a typical self- 
gravitating collisionless galactic disc. Indeed, in the context 
of the upcoming GAIA mission, this externally induced sec¬ 
ular evolution is thought to be a compelling approach to 
describe the radial migration of stars and its impact on 
the observed metallicity gradients (Sellwood & Binney 2002; 
Roskar et al. 2008; Schonrich & Binney 2009; Solway et al. 
2012; Minchev et al. 2013). It may also be applied to de¬ 
scribe the secular diffusion of accretion streams within the 
Galactic halo. Finally, an extension of the formalism of Sec¬ 
tion 3 to discs with a finite thickness might allow us to un¬ 
derstand the process of disc thickening. 
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APPENDIX A: STATISTICAL APPROACH VIA 
HAMILTON’S EQUATION 


We now derive the statistical expression (36) of the sec¬ 
ular diffusion coefficients using a different method based 
on Hamilton’s equations and inspired from Binney & Lacey 
(1988). We will indeed quantify the temporal rate of 
change of the actions, represented by J. The main dif¬ 
ference between the following calculation and that made 
in Binney & Lacey (1988) is that we take explicitly into 
account the self-gravity of the system which leads to the 
appearance of a self-perturbing potential tp” triggered by 
. Starting from the Hamiltonian introduced in equa¬ 
tion (1), and using the Fourier development in angles as 
in equation (8), Hamilton’s equations, 6 = dH/dJ and 
j = —dH/dO, take the form 

= n + , 

< . ^ ^ (Al) 

J = -i^me [tpln +^m] • 

V m 

As we aim to describe the wandering in action-space of the 
particles, we introduce a limited development of the change 
in actions and angles of the form 


e{t) = 00 -I- Ht-I- A0(t), 
J{t) = Jo + A J(t). 


(A2) 


In order to solve this system of coupled differential equa¬ 
tions, we will proceed step-by-step, by including gradually 
the perturbative terms in ip‘^+'ip‘. First of all, one must note 
that the unperturbed orbits follow the straight-line trajec¬ 
tories (0, J) = (00 + Ht, Jo). Then, the first-order term in 
action A J is given by 

AJ(T)= [ dtj{t), (A3) 

Jo 

where J is given by Hamilton’s equations (Al), where all the 
occurences of 0{t) and J{t) are replaced by the expressions 
obtained for the unperturbed orbits. After a time T, the 
shift in action at first order is therefore given by 

AJ{T) = -i^mf dt [V>^(Jo,t)-|-V’^(Jo,t)]e'’^'^®°’''“‘l 

rr,. Jo 


We introduce the operation of angle-average on the initial 
phase 00 as 

= (A4) 

In order to characterize the wandering in action-space, one 
has to study the behavior of the square of the perturbation 
AJ. From Binney & Lacey (1988), we know the relation be¬ 
tween the wandering AJ in action-space and the diffusion 
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coefficient appearing in the Fokker-Planck equation, which 
is given by 


1 


= (T), 


(AS) 


2y L “Jj So' 

where the diffusion equation has been written under the 
compact form 


dt ^ dJi 




(A6) 


Taking an average over the initial phases Oq, using the 
fact that = [V’m**]* and projecting the result on the 

biorthogonal potential basis , one can write 

{A J, A = E E (Jo) X (A7) 

m p,q 

f dti f dt2[ap{ti) + bp{ti)] + . 

Jo Jo 

In order to have an expression which only depends on the 
exterior potential tp^, we use the convolution relation (20), 
written as an amplification relation, to obtain 

{A Ji A = E E E ^ 

m p,q k,l 

f dti f dt 2 fdn 6fc(ri)x (A8) 


/o ^0 
r^2 


(ti—£ 2 ) 


JdT2 (t2-r2)6*(r2)e*'^ “ 

This four dimensional integral is transformed using the 
change of variables 

Ml = ti—ri ; M2 = t2—T2 ; Ml = ri+r2 ; V 2 = T 1 —T 2 ■ (A9) 

It is straightforward to check that the Jacobian of this trans¬ 
formation is 2, so that equation (A8) becomes 

{AJi AJj}g^= E E E ^ 


m p,q k,l 


- / dMi (Mi)e 


im-Cl u\ 


(AlO) 


-im-Cl U 2 


/O 


fT-ui 

dM2 X 
-{T-U2) 


dM 2 (M 2 )e' 

dMi 


f‘2T — ui—U2 — \v2+ui—U2\ 


J\v2\ 

We introduce as in equation (27), the operation of ensemble 
average over many realizations denoted with (.). As in equa¬ 
tion (28), we assume that the exterior perturbing potential 
is a stationary random process. One can then perform an 
ensemble average of the expression (AlO), while assuming, 
as in equation (30), that the response matrix coefficients can 
be taken out of the ensemble average operation. We obtain 

({AJi AJjjg^) = ^^^mirnj %pff [Jo)%pl^* {Jo) 

m p,q k,l 

1 

2 


dMi [I-M]pfc(Mi)e' 


im-Cl u\ 


dM 2 I [I—M] I (M 2 )e ^ 


/ dM2 Cfci(M2) X (All) 

J-(T-U2) 

(2T— Ml— M 2 —|M 2 -|-Mi— M2I —|M 2 |) . 


The next important step of the calculation is to compare T 
with the various autocorrelation times of the system. The 
first one is T^rr describing the typical autocorrelation time 
of the realizations of the external perturbations. Two val¬ 
ues of the potential perturbations separated by a time larger 
than can be considered as independent. The second au¬ 
tocorrelation timescale is Tff,,.,., which describes the typical 
autocorrelation time of the response matrix M and could be 
called the look-back time. From the expression used in equa¬ 
tion (A8), one can note that the values of the self-response 
coefficients are obtained via a non-Markovian mechanism, 
where the past values are amplified thanks to the response 
matrix. However, the self-gravitating system can not have 
an infinite memory, so that only the sufficiently recent past 
values should play a role in this amplification. As a con¬ 
sequence, during the amplification process, only the past 
behavior for a time interval of the order of T^rr is relevant 
and amplified, so that represents the depth with which 
the self-response mechanism can probe past values. We fi¬ 
nally suppose that the time T for which the wandering in 
phase-space is studied satisfies the comparison relations 

; T->TZr- (A12) 

As a consequence, the integration boundaries appearing 
in (All) become 

({A J AJ,}^^) p,^Z\Jo) i’if^iJo) 

m p,q k,l 

1 rT^ 

\ / -^corr — 

dMi[I-M];,i(Mi)e*"‘'^“i X 

^ “dM2 [[I-M]-i]‘(M2)e-^"‘'““= X 

/■Tctrr , _ 

/ dM2 Cfci(M2) e*’” X (A13) 

J-T't 

{2T -Ul-U2-\V2+Ul-U2\-\V2\) ■ 

Thanks to the assumptions (A12), one can see that the last 
term of equation (All) can be approximated by 2T. The re¬ 
maining integrations can then be seen as truncated temporal 
Fourier transforms, so that equation (All) becomes 

({AJAJ,}J = 

T^'^'^mimjpj)^{Jo)'tp)^*{Jo)'S (A14) 

m p,q k,l 

[l-M]-'(m.n)[[l-M];/(m.n)]*CfcKm-n). 

The last step of the simplification is to recall that equa¬ 
tion (19) guarantees that M*=M*, so that using equa¬ 
tion (A5), we finally obtain the expression of the diffusion 
coefficients from equation (A6) which read 

Dij(Jo) =^^^mirrij-tpfi\jo)'ipif*(Jo) x 

m p,q 

C • [I-M]"^l (m-fl) . (A15) 

L ipq 

With this approach, we recover the same diffusion coeffi¬ 
cients as the ones obtained in equation (36) via the quasi- 
linear approach presented in the main text. 
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APPENDIX B: WKB RESPONSE MATRIX 

We estimate the value of the diagonal response matrix coeffi¬ 
cients introduced in equation (69) within the WKB approx¬ 
imation. Using the definition of h{Rg) from equation (70) 
and the fact that is an increasing function of Rg, the 
integration on in equation (69) takes the form 


dRg h{Rg) 


y^27r(cr/-y2)2 


exp 


{R,-Ror 

2{a/V2)^ _ 


(Bl) 


One should note that in this expression, we have a Gaussian 
of spread a/y/2 in Rg, correctly normalized in order to have 
an integral over Rg equal to 1. Assuming that this Gaussian 
is sufficiently peaked, we may replace it by SniRg — Ro), so 
that the integral on J,/, can be dropped. For a Schwarzschild 
distribution function as in equation (45), we then obtain 


M 




(27r)M" 


dJrf, 


dRg 




Rq 


V ^ X 

p LO — nirK — k^Q, 


K , d 
af dJ^, 


ln(^) JdJre ^jl^{Hu^{kr)) 


d 

’dJA 


K 


id Jr JrC ^ J^riHkJkr)) } . (B2) 


In order to perform the integration on Jr, the first step is to 
notice that the only dependence on Jr in the Bessel terms 
is in Hmjykr), through A — \/2JS/k . Therefore, we will 
use the two integration formula (see formula (6.615.1) from 
Gradshteyn & Ryzhik (2007)) 


dJre-’^^" Jl^ipy/Tr) = 




U + oo 

dJrJre-°^^-^jl^{p-FTr) = 

/3^., r/3^1, 


— —- \-\-\-\mr 

2a 


+ 2Q,^l™r|-|-l 


(B3) 


where a>0, /3>0, and nir G E. First, let’s write out explic¬ 
itly the dependence of Hkjykr) with Jr- From equations (43) 
and (65), we find that 


Hk,,{kr) = VXP , (B4) 


where /3 is defined as 


P 


2 

K 


}c2\U2 


2n 1^ 

KRg 



(B5) 


This approximate expression has been obtained using the 
same approximation as in equation (66). We also introduce 
the notation 


We are now able to compute the integrals on Jr from equa¬ 
tion (B2), to obtain 


M 


k,p,kr,Ro\,\k,j,,kr,Ro\ 


{27TfA^ 


dJ^ 


dRg 


UE 


Rq 


E- 




K , d 
-mr^ + k^-— 
ai Ojj, 


(jj — mrK — k,j,n 

/HE 


In - 


\K(T:^ 


(B7) 


‘^dJ. 


^^[P + \mr\-x)Rmr[x]+XR\rr,r\ + l[x\^ ■ 


In order to simplify this expression, we recall that we 
have the property X-mrix)—Rmrix)■ Therefore, in equa¬ 
tion (B7), we have to study four types of sums on rtir, which 
may be simplified as 


V ' TTlr Xm^ (xj 2 ^ , X-gir \x] 

Uj' — rrirK K ^ 1— IcU'/W-rKl^ ’ 

mrSZ mr = l L / ' J 

Rrar [x] ^ _|_ A ^ [x] 

Uj'—nirK LO' LO' 1 — [rTlrK./cj'l^ ’ 

777^’ = ! *■ ' ■* 

E |w.r|Tmr[x] _ 2 Ulr [x] 

ui' — nirR 1 —’ 

771^’ = ! *■ ■’ 

^|m,.| + l[x] _ X\ [x] ^ 2 Xmr + l[x\ 
^ Lo'— rrirR oj' oj' ^ 1 — [mr^/ca'l^ 


(B8) 

(B9) 


where we use uj' =a; —We define the dimensionless 
parameter s as 


uj — kthQ, 

s =- - — 

K 


(BIO) 


We also introduce the reduction factor J^(s,x) (Kalnajs 
1965; Lin & Shu 1966) and similar functions Q{s, x) , R{s, x) 
and X{s, x) defined as 


R{s,x) = 2(l-s^)— 


Xmr [X] 

^ r^A-[s/inrY 


g{s,x) = 2{l-sY- 


X 


1 [x] , 1 Rmr [xl 


2 s 








(Bll) 


T(s,x) = 2(l-s") 


2-,e 


1^1 [X] I 1 ^ 'Rm,r + \-[x\ 


2 s 






Moreover, we notice that we can use the simplification 
d/dJ^K/a'Y\ Or/X =—d/dJ^aX r\^ and that thanks to 
equation (42), one can also explicitly compute 


dJ^, 

RlX 

RoX 

dRg 

Ro 

Ro 2^2 


(B12) 




■ 2fl ■ 

2 

KRg 



(B6) 


Finally, using the expression of the amplitude of the basis 
potentials from equation (63), we obtain a detailled expres- 
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sion of the matrix coefficients as 


AT z' \ ‘^'^GTl\kr\ 


Hs,x) + kl 


K dJtj, [“* \Ka^ 


In 


ns 


Sis,x) 


(B13) 


_j_ 


[{'^-x)Q{s,x)+'H{s,x)+X^{s,x)] 


where one must remember that within the WKB approxi¬ 
mation, the response matrix is diagonal. For a tepid disc, we 
may neglect some of the terms appearing in equation (B13). 
A tepid disc corresponds to a disc where the orbits pos¬ 
sess a small radial energy, so that all the orbits are close to 
circular orbits. It also implies that \dFo/dJr\ S> 

For a Schwarzschild distribution function, the typical spread 
in Jr is of the order of a^jn, so that we may consider 
equation (B13) as a limited development in and 

d// k\. Therefore, for a tepid disc the diagonal co¬ 
efficients of the response matrix finally take the form given 
in equation (74). 


APPENDIX C: AUTOCORRELATION 
DIAGONALIZATION 

Let us now show how the hypothesis of spatially quasi- 
stationarity of the external perturbations introduced in 
equation (98) leads to a diagonalization of the autocorre¬ 
lation with respect to the radial frequencies kr as shown in 
equation (99). In order to shorten the notations, we do not 
write anymore the dependence with respect to the azimuthal 
number m,/,, and the exterior perturbation will be noted as 
i/) = 'i/j®. As a consequence, the assumption of temporal and 
quasi-spatial stationarity from equation (98) takes the form 

{tp[Rl,tl]tp*[R2,t2]) = c[tl—f2, Rl—R 2 , (Ri+R 2)/2] . (Cl) 

Equation (94) for the diffusion coefficients requires us to 
study the term [i?g, wi] ^^*2 [Rg,CJ 2 ]). Using the defini¬ 
tion of the temporal Fourier transform from equation (11) 
and the local radial Fourier transform from equation (84), 
we may rewrite it as 

(^fcl[^9,^l] = (C2) 

-^-jJdtidt 2 dRidR 2 ti] '!/’*[R 2 , ^ 2 ]} x 

g[Rg-Ri] g[Rg-R 2 ] , 

where g\R] is defined as 

fif[R] = exp [-R^/(2cr^)] . (C3) 

We now use the assumption from equation (Cl) relative to 
the radial dependences of the perturbation autocorrelation, 
and the change of variables 

\ut=tl+t2 ; Vt=tl—t2, 

s , (C4) 

Mr = 2 (^ 1 +^ 2 ) ; Mr = Rl — R 2 ■ 


This transformation is of determinant 2, so that equa¬ 
tion (C2) becomes 


j dut dvt dUr dVr 6* 


> . 0^1 -i-ojn 




-Vr — ^ 


(C5) 


5r[Rg—Mr —Mr/2] g [Rg — Mr-|-Mr/2] C [Mt, Mr, Mr] . 


The integration on ut is straightforward and is equal to 
27r(5D((tMi— cJ 2)/2). The integration on vt is then direct and 
gives C [cMi,Mr,Mr]. Finally, we note that the product of the 
two Gaussians in equation (C5) can be rewritten in order to 
disentangle the dependences on Mr and Mr to read 

5l[Rg-Mr-Mr/2]g(Rg-Mr-|-Mr/2] =g(v^(Rg-Mr)]s]Mr/v^] , (C6) 

where the presence of ^/2 comes from the definition of the g 
function introduced in equation (C3). One can then rewrite 
equation (C5) as 


V'fc2[-R9,tM2]) = (C7) 

^e'«»('='-'=')5D(cMl-CM2)y’dMre-'^’'’'3[Mr/V2] X 


IdUr C [iOl ,Vr,Ur]g[V2(Rg-Ur)]e 




As we have assumed that the function Ur C [wi, Mr, Mr] is 
a slowly varying function, we may take it out of the integra¬ 
tion on Mr and evaluate it as C [wi, Mr, Rg]. The remaining 
integration on Mr can then be computed and reads 


(dMr(7[v^(Rg-Mr)]e-*''''-'''>“’' 

^ — iRg(k^ — k‘^) 


= y/nae 


exp 


(C 8 ) 

{kl-klf 


4/(72 






where we replaced the Gaussian in fc/ —fcr by a Dirac delta, 
while paying a careful attention to the correct normalization. 
As a consequence, equation (C7) becomes 


((/)i,l[Rg,Wl] V'fc*2[R9,CM2]) = (5d(wi-W 2) 5D(fc/-fer) X (C 9 ) 


JdVr e [wi, Mr, Rg] g[Vr/V2] . 


Because of the definition from equation (84), the presence 
of the factor Xjypl corresponds to the change a —>■ ypl a so 
that the remaining integral on Vr may be interpreted as a 
local radial Fourier transform centered around the position 
Mr = 0. Therefore, we straightforwardly obtain the diagonal¬ 
ized expression introduced in equation (99). 



